Mixing times in evolutionary game dynamics 
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Without mutation and migration, evolutionary dynamics ultimately leads to the extinction of all 
but one species. Such fixation processes are well understood and can be characterized analytically 
with methods from statistical physics. However, many biological arguments focus on stationary 
distributions in a mutation-selection equilibrium. Here, we address the equilibration time required 
to reach stationarity in the presence of mutation, this is known as the mixing time in the theory 
of Markov processes. We show that mixing times in evolutionary games have the opposite be- 
haviour from fixation times when the intensity of selection increases: In coordination games with 
bistabilities, the fixation time decreases, but the mixing time increases. In coexistence games with 
metastable states, the fixation time increases, but the mixing time decreases. Our results are based 
on simulations and the WKB approximation of the master equation. 

PACS numbers: 02.50.Ga Markov processes, 05.10.Gg Stochastic analysis methods (Fokker-Planck, Langcvin, 
etc.), 89.75.-k (Complex systems), 03.65.Sq (Semiclassical theories and applications) 



How long does it take for a stochastic many-particle 
system to reach its stationary distribution? This question 
goes beyond traditional equilibrium statistical physics 
and requires a theory for non-equilibrium systems. Sig- 
nificant progress has been made over the last decades, 
but developing a more complete theory is still very much 
work in progress. Many non-equilibrium systems lack an 
energy or Lyapunov function, any theoretical analysis has 
to start from the microscopic dynamics itself. Such ap- 
proaches have been applied successfully to off-equilibrium 
phenomena in physics [1] , but also to a number of appli- 
cations in adjacent disciplines, including in epidemiol- 
ogy, biological transport, pattern formation, agent-based 
models in economics and of social phenomena [2]-[5] . 

For stochastic processes with absorbing states, our 
opening question can be answered - at least to some ex- 
tent. Absorbing states are those in which the system gets 
'trapped', so that a full dynamic arrest occurs. Systems 
with absorbing states exhibit new types of phase tran- 
sitions, universality classes and complexity, previously 
unknown in physics [6l [7]. They are relevant in social 
systems, where an absorbing state may correspond to 
a uniform consensus, and in evolutionary biology where 
they describe fixation of a trait. Stochasticity can also 
drive individual phenotypes to extinction in evolutionary 
game dynamics. In the absence of mutation, a given phe- 
notype is never re-introduced once it has become elim- 
inated from the population. Answering the question of 
equilibration times for this type of systems then amounts 
to calculating the time to fixation [51 [5] ■ 

The purpose of our work is to develop a similar ap- 
proach for evolutionary systems with mutation. In such 
systems there are no absorbing states and thus no fixa- 
tion. Still, they reach a stationary distribution at asymp- 
totic times. In order to characterize the approach to sta- 



tionarity we consider what is referred to as the mixing 
time in the theory of Markov process [TO]. This is the 
time needed for the probability distribution over states 
to approach its stationary shape up to some specified 
small distance s. Mixing times have been considered in 
the context of Markov Chain Monte Carlo methods [TU] , 
and recently in game dynamical learning [TTJ, but, to 
our knowledge, they have not been discussed for evolu- 
tionary processes. We here introduce the basic concepts, 
analyze mixing times in 2 x 2 evolutionary games and 
show how methods from quasi-classical physics can be 
used to obtain analytical approximations. Our analysis 
is based on computer simulations and analytical approx- 
imations using the Wentzel-Kramers-Brillouin (WKB) 
method [12j [24]. While we focus on specific instances 
of evolutionary dynamics, we expect that these tools can 
be applied to describe the non-equilibrium dynamics of 
a large class of individual-based models. 

We consider a well-mixed population of iV individuals 
of type A or B. The state of the system is determined 
by the number n £ {0, . . . , N} of individuals of type A. 
The fitness of individuals of the two types interacting in 
an evolutionary 2x2 game is given by [T3] 

n B (n) =^c+^^d, (1) 

if the system is in state n. The parameters a, 6, c and d 
specify the underlying game. We study the evolutionary 
dynamics defined by the birth-death process with rates 

m 

T~ = ^[l- J 9AII(n)]=^ + £* (2) 



2 



where AII(n) = 11^ (n) — Ils(n), and where T+ is the 
rate at which a individual of type B is replaced by an 
individual of type A in state n, T~ is the rate of the 
opposite event. The parameter [i S [0, 1] represents the 
mutation rate, (3 > is the intensity of selection. The 
probability Pt(n) of finding the system in state n at time 
t then obeys the master equation 

P t (n) = T+_ 1 P t (n-l)+T- +1 P t (n+l)-(T++T-)P t (n). 

(3) 

We will denote the stationary distribution at asymp- 
totic times by ip*(n) — lim t _>. 00 P t (n) . Following 
[TU] . the mixing time t m i x (s) is defined as t m i x (e) = 
min{t : d(t) < e} where the variational distance d(t) = 
\ T n \Pt{n) — ili*(n)\ measures the distance of the prob- 
ability distribution Pt(-) from the stationary distribution 

r(-) m- 

In order to characterize the behaviour of mixing times 
in different scenarios, it is useful to first consider the 
limit of N — > oo. In this case, the fraction of individuals 
of type A, x = n/N , fulfills the deterministic replicator- 
mutator equation (RME) 

x = /3x(l - x)[(a-c)x - (d-b)(l - x)](l - //) + ^M 4 ) 

The number, position and stability of the fixed points of 
Eq. Q generally depend on the parameters (3 and u, as 
well as on the underlying game [T31 [IS] . 

We first study a symmetric coexistence game, defined 
by a = 1, b = 2, c = 2, d = 1. In this scenario, Eq. Q 
has one stable fixed point at x* = 1/2, and no other 
fixed points. The left panel of Fig. [TJsliows the stationary 
distributions of the resulting Markov chain for different 
intensities of selection (3. With increasing /3, the distri- 
bution concentrates on the area around the deterministic 
fixed point. The corresponding mixing times are shown 
in the right-hand panel, starting the dynamics in a sin- 
gle state, Po(n) = <5„,„ , for varying no- Increasing the 
intensity of selection f3 reduces the mixing time, because 
mixing is governed by the deterministic flux, which in- 
creases with (3. For no ^ the mixing time is limited 
by this term. Note that the fixation time in this game 
increases exponentially with N, as motion against the 
deterministic flux is required [3]. 

Next, we address a symmetric coordination game with 
parameters a = 2, b = 1, c = 1, d = 2. In this case, 
Eq. Q has an unstable fixed point at x\ = \ and two 
stable fixed points x-i and X3 near x = and x = 1 respec- 
tively. Fig. [2] shows the bimodal stationary distributions 
and mixing times for this game. With increasing /?, the 
distribution becomes sharply peaked around the two sta- 
ble fixed points. The mixing time is then very sensitive to 
the initial condition. When the system is started near to 
either of the stable fixed points, the system will quickly 
tend to a quasi-stationary distribution (QSD) around one 
of the fixed points. The probability will then slowly leak 
over to the other side on a time scale exponentially slow 
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FIG. 1: (Colour online) Symmetric co-existence game in a 
population of N — 100 individuals. Left: Stationary distri- 
butions of the stochastic dynamics (Eq. (|3|). Right: Mixing 
time (e = 1/4) when the stochastic process is started from a 
Dirac distribution at Po{n) = S(n — no). The mutation rate 
is u — 1/101, leading to a uniform stationary distribution at 
13 = 0. 

in N. Thus, the mixing times increase with increasing 
intensity of selection (3, whereas fixation times would de- 
crease with (3 [5] . We can exploit this separation of time 
scales to calculate the mixing time analytically, similar 
to the problem of calculating the mean switching time 
between quasi-stationary states. 

Let us assume that we start to the left of the unstable 
fixed point. The time-dependent ansatz we use for the 
probability distribution is 

V [n) ~ \ ^*(n)(l - e~ Et ) n>m' [) 

where — E < is the eigenvalue of the slowest decaying 
mode of the problem, and where n\ corresponds to the 
central unstable fixed point. This ansatz is valid on an 
exponentially long time scale in N. Explicitly calculating 
the variational distance between the ansatz of Eq. ([5]) and 
the stationary distribution, we find i mix (e) = — _E~Mn[2e] 
such that the problem reduces to finding the eigenvalue 
— E. Based on Eq. the current through the cen- 
tral fixed point, is given by J(t) — ^ Yl n <m ip lcllk ( n ) — 
— (E/2)e~ Et . Thus we can find E from the initial cur- 
rent, J(Q) = —E/2. Calculating the mixing time then 
reduces to the well studied problem of determining the 
escape current in a bistable potential |16) . 

As this is a one-step process, exact expressions exist 
for the mean first passage times |17j . One avenue would 
be to derive the large ./V asymptotics for these O HH [19] . 
We do not follow this approach here, instead we calcu- 
late the initial current through the unstable fixed point 
based on the celebrated WKB approximation. This has 
two advantages: first, this method is valid for a much 
wider range of problems, such as those with multiple 
jumps [201 HI]: or of higher dimensions |2]Jj. Secondly, 
the stationary distribution ip*(n) is calculated as a by- 
product. Our approach also complements recent stud- 



3 



0.06 p. 
0.05 - 



,8=0.75 




(i.i)o L 



20 40 60 80 




20 40 60 80 100 
Initial condition no 



FIG. 2: (Colour online) Symmetric coordination game in a 
population of N = 100 individuals. Left: Stationary distribu- 
tion of the master equation pj. Right: Mixing time (e = 1/4) 
as a function the position of the Dirac distribution from which 
the dynamics is started. The mutation rate is /i = 1/101, 
leading to a flat stationary distribution for /3 = 0. 



and the relaxation solution S%(x) = \n[H p (x,0)], where 
H p = dH/dp. In our setup the activation solution de- 
scribes the behaviour of the QSD to the left of x% = 1/2. 
The relaxation mode, describing deterministic motion to 
the right of x\, will play no significant role in our further 
analysis. 

In order to complete the calculation two main tasks 
remain: (a) the activation solution ip WKB defined by Eqs. 
( 6|8|9 1 needs to be normalized, and (b) the QSD needs 
to be connected to the initial current J(0) through X\. 
These tasks can be addressed by performing a Kramers- 
Moyal expansion of the master equation ^ around the 
unstable fixed point x\. Writing f±(x) — Nw±(x)ip(x) 
we find 



r=±l 



2A 2 



(10) 



ies which have successfully introduced these methods to 
evolutionary game theory by calculating fixation times 
in evolutionary games without mutation |23j . Before we 
describe the main steps of the calculation, it is useful to 
introduce x = n/N and to expand the transition rates, 
Eq. ([2]), into powers of 1/N in leading and sub-leading or- 
der. Specifically, we write (Nx) = w± (x) + u± (x) /N 
[2"T] . We use the WKB ansatz, 



f KB (i) = cxp[-iVS(a:) - S x (x)] 



(6) 



where both S{x) and S\(x) are assumed to be of order 
./V . It is important to note the difference between ansatz 
([5]) and (JoJ) . Ansatz ^ is valid on exponentially long time 
scales and about the asymptotes of the distribution and 
it takes account of the back current through the central 
fixed point. Ansatz ^ is valid everywhere, but only on 
short time scales. It can therefore be used to calculate 
the initial current. We proceed by inserting Eq. Mm into 
Eq. (J3J , assuming (quasi-) stationarity (dtip WK = 0) 
and expanding the resulting equation into powers of A? --1 
[101 HD 123 122- In lowest order, we find a Hamilton- 
Jacobi equation, 

H(x,p) = w + {x){e p - 1) + w-(x)(e- p - 1) = 0, (7) 

where p = d x S. This constitutes an equation for S(x) 
and it has two solutions: (i) the activation solution 



S(x) 



d£ln 



«>+(0 



(8) 



and (ii) the so-called relaxation solution S(x) = 0. In 
next order, we find the activation solution 



Si(x) 



«-(0 _ 



1 



hi[w + (x)w_(x)], 
(9) 



where the term in the square bracket is identified as the 
divergence-free probability current J(0). Further alge- 
braic manipulations then lead to [20] . 

H pp (x,0)d x iP WKB (x) 



2N 



(11) 



J = Tjj WKU (x)H px (x 1 0)(x-x 1 ) 
where H pp — d 2 H/dp 2 = 

d 2 H I dxdp = W—w',/w+ — w+w'_/w-. Re-arranging Eq. 
(Ill one has 



w_ and H px — 



6 

H px £ 



crfc(y), where y = 



x — X\ 



\J Hpp/ (N H px ) 
(12) 

The final step then consists in matching the asymp- 
tote, tp WKB (y) = J^e y2 , valid for y < -1, with the 
relaxation-mode solution ip WKB (x) — Ae~ NS ^~ Sl ^ 
(with S(x) and Si(x) given by Eqs. ^ and ^). The 
normalisation constant A is obtained from a Gaussian 
approximation of the relaxation solution about the fixed 
point X2 [22 . 

Carrying out this procedure the initial current is found 
to be exponential in N: 



J(0) 



Z^t/S»( X2 )\S»{ Xi )\ 



xe 



N[S(x 2 )-S(x 1 )]+S 1 (x 2 )-S{x 1 ) 



Finally the mixing time is 



t mix (e) = ln(2 £ )/(2J(0)). 



(13) 



(14) 



Figure [3] shows the mixing times calculated via the 
WKB method along with direct computation from a nu- 
merical integration of the master equation. Agreement is 
generally very good, except for small values of N when 
the introduction of continuous variable x as well as the 
expansion in powers of A -1 become inaccurate. The 
slight offset between the two sets of results is due to the 
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Population size N 

FIG. 3: (Colour online) Mixing time (e = 1/4) for the sym- 
metric coordination game as a function of the population size 
AT. Lines are calcuated from Eq. ( |14[ ), dots are from integra- 
tion of the master equation. 



error introduced in the calculation by assuming Gaus- 
sian approximation made when normalizing the QSD 
-0 WKB (-). Better agreement could be obtained by nor- 
malizing the distribution exactly. However, this requires 
a numerical approach, whereas our final result is more 
explicit. 

While the WKB approach can successfully be em- 
ployed to obtain mixing times, there are limitations to 
this method. One potential problem is the divergence of 
the WKB solution i/> WKB (cc) at the boundaries of the sys- 
tem. This does not affect the outcome of our calculation 
as long as the stable fixed points of the RME are not too 
close to the boundaries of phase space, but it does limit 
the range of j3 and u for which it is valid. We stress that 
the ansatz ^ still applies, but the eigenvalue E needs to 
be calculated via a different approach. The methods we 
have presented lend themselves to generalisation. For ex- 
ample the assumption of symmetry of the problem can be 
relaxed, the unstable fixed point need not be at x = 1/2. 

In summary, we have introduced the concept of mix- 
ing times for evolutionary dynamics with mutation. As 
intensity of selection is increased, the mixing times in 
co-existence games decrease. In coordination games, one 
observes the opposite trend. In both cases the behaviour 
of mixing times is opposite to that of fixation times in 
the corresponding systems without mutation. The con- 
cept of mixing times may often be more appropriate for 
many biological systems than the computation of fixa- 
tion times, in particular when effects of mutation or im- 
migration cannot be ignored [25] . As shown in our work, 
tools from theoretical physics can be used to successfully 
estimate mixing times based on semi-analytical consider- 
ations. We expect this to be useful not only for biological 
systems, but also for models of social dynamics and other 
interacting many-particle processes. 
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